library(data.table)
library(ggplot2)

base <- "/data/workspace/sanction/aamas15-code/output/runs-stability/runs"
dirs <- c(5, 10, 20, 30, 50, 100)
contractF <- "contract.csv"
prosumerF <- "prosumer.csv"
providerF <- "provider.csv"
regulatorF <- "regulator.csv"
cycles <- 1000

##
## PRE-PROCESSING
##
for(dir in dirs) {
  
  mergeNV <- data.table(cycle=1:cycles-1)
  setkey(mergeNV, cycle)

  mergeNC <- data.table(cycle=1:cycles-1)
  setkey(mergeNC, cycle)

  mergePV <- data.table(cycle=1:cycles-1)
  setkey(mergePV, cycle)

  mergeAPUN <- data.table(cycle=1:cycles-1)
  setkey(mergeAPUN, cycle)

  mergeAINFO <- data.table(cycle=1:cycles-1)
  setkey(mergeAINFO, cycle)

  mergeAPINFO <- data.table(cycle=1:cycles-1)
  setkey(mergeAPINFO, cycle)
  
  mergeRPUN <- data.table(cycle=1:cycles-1)
  setkey(mergeRPUN, cycle)
  
  mergeRINFO <- data.table(cycle=1:cycles-1)
  setkey(mergeRINFO, cycle)
  
  mergeRPINFO <- data.table(cycle=1:cycles-1)
  setkey(mergeRPINFO, cycle)
  
  mergeADEN <- data.table(cycle=1:cycles-1)
  setkey(mergeADEN, cycle)
  
  mergeAREP <- data.table(cycle=1:cycles-1)
  setkey(mergeAREP, cycle)
  
  mergeRREP <- data.table(cycle=1:cycles-1)
  setkey(mergeRREP, cycle)
  
  mergeRSUS <- data.table(cycle=1:cycles-1)
  setkey(mergeRSUS, cycle)

  mergeQS <- data.table(cycle=1:cycles-1)
  setkey(mergeQS, cycle)

  mergeQB <- data.table(cycle=1:cycles-1)
  setkey(mergeQB, cycle)

  replicas <- dir
  for(replica in 0:(replicas-1)) {
    ##
    ## CONTRACT
    ##
    filename <- paste0(base,dir,"/",replica,"/",contractF)
    contractD <- data.table(read.csv(filename, header=TRUE, sep=";",
                                   quote="\""))
    setkey(contractD, cycle)
  
    aux <- contractD[, numV:=nrow(.SD[which(sellerComplied == "false")]), by=cycle]
    aux <- contractD[, numC:=nrow(.SD[which(sellerComplied == "true")]), by=cycle]
    aux <- contractD[, prop:=nrow(.SD[which(sellerComplied == "false")]) / .N, by=cycle]
    
    meanNV <- aux[, mean(numV), by=cycle]
    setnames(meanNV, c("cycle","V1"),c("cycle",paste0("V",replica)))
    setkey(meanNV, cycle)
    mergeNV <- merge(mergeNV, meanNV, all = FALSE)
    
    meanNC <- aux[, mean(numC), by=cycle]
    setnames(meanNC, c("cycle","V1"),c("cycle",paste0("V",replica)))
    setkey(meanNC, cycle)
    mergeNC <- merge(mergeNC, meanNC, all = FALSE)
    
    meanPV <- aux[, mean(prop), by=cycle]
    setnames(meanPV, c("cycle","V1"),c("cycle",paste0("V",replica)))
    setkey(meanPV, cycle)
    mergePV <- merge(mergePV, meanPV, all = FALSE)
    
    ##
    ## PROSUMER
    ##
    filename <- paste0(base,dir,"/",replica,"/",prosumerF)
    prosumerD <- data.table(read.csv(filename, header=TRUE, sep=";",
                                     quote="\""))
    setkey(prosumerD, cycle)
    
    aux <- prosumerD[, apun:=sum(.SD$a_pun_sanction), by=cycle]
    aux <- prosumerD[, apinfo:=sum(.SD$a_pun.info_sanction), by=cycle]
    aux <- prosumerD[, ainfo:=sum(.SD$a_info_sanction), by=cycle]
    aux <- prosumerD[, rpun:=sum(.SD$r_pun_sanction), by=cycle]
    aux <- prosumerD[, rpinfo:=sum(.SD$r_pun.info_sanction), by=cycle]
    aux <- prosumerD[, rinfo:=sum(.SD$r_info_sanction), by=cycle]
        
    meanAPUN <- aux[, mean(apun), by=cycle]
    setnames(meanAPUN, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeAPUN <- merge(mergeAPUN, meanAPUN, all = FALSE)
    
    meanAPINFO <- aux[, mean(apinfo), by=cycle]
    setnames(meanAPINFO, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeAPINFO <- merge(mergeAPINFO, meanAPINFO, all = FALSE)
  
    meanAINFO <- aux[, mean(ainfo), by=cycle]
    setnames(meanAINFO, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeAINFO <- merge(mergeAINFO, meanAINFO, all = FALSE)
    
    meanRPUN <- aux[, mean(rpun), by=cycle]
    setnames(meanRPUN, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeRPUN <- merge(mergeRPUN, meanRPUN, all = FALSE)
    
    meanRPINFO <- aux[, mean(rpinfo), by=cycle]
    setnames(meanRPINFO, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeRPINFO <- merge(mergeRPINFO, meanRPINFO, all = FALSE)
    
    meanRINFO <- aux[, mean(rinfo), by=cycle]
    setnames(meanRINFO, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeRINFO <- merge(mergeRINFO, meanRINFO, all = FALSE)
    
    aux <- prosumerD[, aden:=sum(.SD$a_den_sanction), by=cycle]
    aux <- prosumerD[, arep:=sum(.SD$a_rep_sanction), by=cycle]
    aux <- prosumerD[, rrep:=sum(.SD$r_rep_sanction), by=cycle]
    aux <- prosumerD[, rsus:=sum(.SD$r_sus_sanction), by=cycle]
      
    meanADEN <- aux[, mean(aden), by=cycle]
    setnames(meanADEN, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeADEN <- merge(mergeADEN, meanADEN, all = FALSE)
    
    meanAREP <- aux[, mean(arep), by=cycle]
    setnames(meanAREP, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeAREP <- merge(mergeAREP, meanAREP, all = FALSE)
      
    meanRREP <- aux[, mean(rrep), by=cycle]
    setnames(meanRREP, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeRREP <- merge(mergeRREP, meanRREP, all = FALSE)
    
    meanRSUS <- aux[, mean(rsus), by=cycle]
    setnames(meanRSUS, c("cycle","V1"),c("cycle",paste0("V",replica)))
    mergeRSUS <- merge(mergeRSUS, meanRSUS, all = FALSE)
  
    ##
    ## PROVIDER
    ##
    filename <- paste0(base,dir,"/",replica,"/",providerF)
    providerD <- data.table(read.csv(filename, header=TRUE, sep=";",
                                     quote="\""))
    setkey(providerD, cycle)
  
    meanQS <- providerD[, mean(quantity_sell), by=cycle]
    setnames(meanQS, c("cycle","V1"),c("cycle",paste0("V",replica)))
    setkey(meanQS, cycle)
    mergeQS <- merge(mergeQS, meanQS, all = FALSE)
  
    meanQB <- providerD[, mean(quantity_buy), by=cycle]
    setnames(meanQB, c("cycle","V1"),c("cycle",paste0("V",replica)))
    setkey(meanQB, cycle)
    mergeQB <- merge(mergeQB, meanQB, all = FALSE)
  }
  
  filename <- paste0(base,dir)
  write.table(mergeNV, file = paste0(filename,"/mergeNV.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeNC, file = paste0(filename,"/mergeNC.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergePV, file = paste0(filename,"/mergePV.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeAPUN, file = paste0(filename,"/mergeAPUN.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeAINFO, file = paste0(filename,"/mergeAINFO.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeAPINFO, file = paste0(filename,"/mergeAPINFO.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeRPUN, file = paste0(filename,"/mergeRPUN.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeRINFO, file = paste0(filename,"/mergeRINFO.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeRPINFO, file = paste0(filename,"/mergeRPINFO.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  
  write.table(mergeADEN, file = paste0(filename,"/mergeADEN.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  write.table(mergeAREP, file = paste0(filename,"/mergeAREP.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE) 
  write.table(mergeRREP, file = paste0(filename,"/mergeRREP.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE) 
  write.table(mergeRSUS, file = paste0(filename,"/mergeRSUS.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
  
  write.table(mergeRSUS, file = paste0(filename,"/mergeQS.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE) 
  write.table(mergeRSUS, file = paste0(filename,"/mergeQB.csv"), sep = ";",
              row.names = FALSE, col.names = TRUE, append = FALSE)
}


##
## ANALYSIS
##
NV <- data.table(cycle=1:cycles-1)
setkey(NV, cycle)

NC <- data.table(cycle=1:cycles-1)
setkey(NC, cycle)

PV <- data.table(cycle=1:cycles-1)
setkey(PV, cycle)

APUN <- data.table(cycle=1:cycles-1)
setkey(APUN, cycle)

AINFO <- data.table(cycle=1:cycles-1)
setkey(AINFO, cycle)

APINFO <- data.table(cycle=1:cycles-1)
setkey(APINFO, cycle)

RPUN <- data.table(cycle=1:cycles-1)
setkey(RPUN, cycle)

RINFO <- data.table(cycle=1:cycles-1)
setkey(RINFO, cycle)

RPINFO <- data.table(cycle=1:cycles-1)
setkey(RPINFO, cycle)

ADEN <- data.table(cycle=1:cycles-1)
setkey(ADEN, cycle)

AREP <- data.table(cycle=1:cycles-1)
setkey(AREP, cycle)

RREP <- data.table(cycle=1:cycles-1)
setkey(RREP, cycle)

RSUS <- data.table(cycle=1:cycles-1)
setkey(RSUS, cycle)

QS <- data.table(cycle=1:cycles-1)
setkey(QS, cycle)

QB <- data.table(cycle=1:cycles-1)
setkey(QB, cycle)

for(dir in dirs){
  replicas <- dir
  filename <- paste0(base,dir)
  
  mergeNV <- data.table(read.csv(paste0(filename,"/mergeNV.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeNC <- data.table(read.csv(paste0(filename,"/mergeNC.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergePV <- data.table(read.csv(paste0(filename,"/mergePV.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeAPUN <- data.table(read.csv(paste0(filename,"/mergeAPUN.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeAINFO <- data.table(read.csv(paste0(filename,"/mergeAINFO.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeAPINFO <- data.table(read.csv(paste0(filename,"/mergeAPINFO.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeRPUN <- data.table(read.csv(paste0(filename,"/mergeRPUN.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeRINFO <- data.table(read.csv(paste0(filename,"/mergeRINFO.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeRPINFO <- data.table(read.csv(paste0(filename,"/mergeRPINFO.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeADEN <- data.table(read.csv(paste0(filename,"/mergeADEN.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeAREP <- data.table(read.csv(paste0(filename,"/mergeAREP.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeRREP <- data.table(read.csv(paste0(filename,"/mergeRREP.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeRSUS <- data.table(read.csv(paste0(filename,"/mergeRSUS.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeQS <- data.table(read.csv(paste0(filename,"/mergeQS.csv"), header=TRUE,
                                 sep=";", quote="\""))
  mergeQB <- data.table(read.csv(paste0(filename,"/mergeQB.csv"), header=TRUE,
                                 sep=";", quote="\""))
  
  mergePV$mean <- apply(mergePV[,colnames(mergePV)[2:replicas+1], with=FALSE], 1, mean)
  mergePV$sd <- apply(mergePV[,colnames(mergePV)[2:replicas+1], with=FALSE], 1, sd)
  PV <- cbind(PV, mergePV$mean, mergePV$sd)
    
  mergeNV$mean <- apply(mergeNV[,colnames(mergeNV)[2:replicas+1], with=FALSE], 1, mean)
  mergeNV$sd <- apply(mergeNV[,colnames(mergeNV)[2:replicas+1], with=FALSE], 1, sd)
  NV <- cbind(NV, mergeNV$mean, mergeNV$sd)
    
  mergeNC$mean <- apply(mergeNC[,colnames(mergeNC)[2:replicas+1], with=FALSE], 1, mean)
  mergeNC$sd <- apply(mergeNC[,colnames(mergeNC)[2:replicas+1], with=FALSE], 1, sd)
  NC <- cbind(NC, mergeNC$mean, mergeNC$sd)
  
  mergeAPUN$mean <- apply(mergeAPUN[,colnames(mergeAPUN)[2:replicas+1], with=FALSE], 1, mean)
  mergeAPUN$sd <- apply(mergeAPUN[,colnames(mergeAPUN)[2:replicas+1], with=FALSE], 1, sd)
  APUN <- cbind(APUN, mergeAPUN$mean, mergeAPUN$sd)
  
  mergeAPINFO$mean <- apply(mergeAPINFO[,colnames(mergeAPINFO)[2:replicas+1], with=FALSE], 1, mean)
  mergeAPINFO$sd <- apply(mergeAPINFO[,colnames(mergeAPINFO)[2:replicas+1], with=FALSE], 1, sd)
  APINFO <- cbind(APINFO, mergeAPINFO$mean, mergeAPINFO$sd)
  
  mergeAINFO$mean <- apply(mergeAINFO[,colnames(mergeAINFO)[2:replicas+1], with=FALSE], 1, mean)
  mergeAINFO$sd <- apply(mergeAINFO[,colnames(mergeAINFO)[2:replicas+1], with=FALSE], 1, sd)
  AINFO <- cbind(AINFO, mergeAINFO$mean, mergeAINFO$sd)
  
  mergeRPUN$mean <- apply(mergeRPUN[,colnames(mergeRPUN)[2:replicas+1], with=FALSE], 1, mean)
  mergeRPUN$sd <- apply(mergeRPUN[,colnames(mergeRPUN)[2:replicas+1], with=FALSE], 1, sd)
  RPUN <- cbind(RPUN, mergeRPUN$mean, mergeRPUN$sd)
  
  mergeRPINFO$mean <- apply(mergeRPINFO[,colnames(mergeRPINFO)[2:replicas+1], with=FALSE], 1, mean)
  mergeRPINFO$sd <- apply(mergeRPINFO[,colnames(mergeRPINFO)[2:replicas+1], with=FALSE], 1, sd)
  RPINFO <- cbind(RPINFO, mergeRPINFO$mean, mergeRPINFO$sd)
  
  mergeRINFO$mean <- apply(mergeRINFO[,colnames(mergeRINFO)[2:replicas+1], with=FALSE], 1, mean)
  mergeRINFO$sd <- apply(mergeRINFO[,colnames(mergeRINFO)[2:replicas+1], with=FALSE], 1, sd)
  RINFO <- cbind(RINFO, mergeRINFO$mean, mergeRINFO$sd)
  
  mergeADEN$mean <- apply(mergeADEN[,colnames(mergeADEN)[2:replicas+1], with=FALSE], 1, mean)
  mergeADEN$sd <- apply(mergeADEN[,colnames(mergeADEN)[2:replicas+1], with=FALSE], 1, sd)
  ADEN <- cbind(ADEN, mergeADEN$mean, mergeADEN$sd)
  
  mergeAREP$mean <- apply(mergeAREP[,colnames(mergeAREP)[2:replicas+1], with=FALSE], 1, mean)
  mergeAREP$sd <- apply(mergeAREP[,colnames(mergeAREP)[2:replicas+1], with=FALSE], 1, sd)
  AREP <- cbind(AREP, mergeAREP$mean, mergeAREP$sd)
  
  mergeRREP$mean <- apply(mergeRREP[,colnames(mergeRREP)[2:replicas+1], with=FALSE], 1, mean)
  mergeRREP$sd <- apply(mergeRREP[,colnames(mergeRREP)[2:replicas+1], with=FALSE], 1, sd)
  RREP <- cbind(RREP, mergeRREP$mean, mergeRREP$sd)
  
  mergeRSUS$mean <- apply(mergeRSUS[,colnames(mergeRSUS)[2:replicas+1], with=FALSE], 1, mean)
  mergeRSUS$sd <- apply(mergeRSUS[,colnames(mergeRSUS)[2:replicas+1], with=FALSE], 1, sd)
  RSUS <- cbind(RSUS, mergeRSUS$mean, mergeRSUS$sd)
  
  mergeQS$mean <- apply(mergeQS[,colnames(mergeQS)[2:replicas+1], with=FALSE], 1, mean)
  mergeQS$sd <- apply(mergeQS[,colnames(mergeQS)[2:replicas+1], with=FALSE], 1, sd)
  QS <- cbind(QS, mergeQS$mean, mergeQS$sd)
  
  mergeQB$mean <- apply(mergeQB[,colnames(mergeQB)[2:replicas+1], with=FALSE], 1, mean)
  mergeQB$sd <- apply(mergeQB[,colnames(mergeQB)[2:replicas+1], with=FALSE], 1, sd)
  QB <- cbind(QB, mergeQB$mean, mergeQB$sd)
}
cNames <- c("cycle",
            "R5m","R5s","R10m","R10s",
            "R20m","R20s","R30m","R30s",
            "R50m","R50s","R100m","R100s")
setnames(NV, cNames)
setnames(NC, cNames)
setnames(PV, cNames)
setnames(APUN, cNames)
setnames(AINFO, cNames)
setnames(APINFO, cNames)
setnames(RPUN, cNames)
setnames(RINFO, cNames)
setnames(RPINFO, cNames)
setnames(ADEN, cNames)
setnames(AREP, cNames)
setnames(RREP, cNames)
setnames(RSUS, cNames)
setnames(QS, cNames)
setnames(QB, cNames)

sd(PV[which(cycle > 800)]$R5m) / mean(PV[which(cycle > 800)]$R5m)
sd(PV[which(cycle > 800)]$R10m) / mean(PV[which(cycle > 800)]$R10m)
sd(PV[which(cycle > 800)]$R20m) / mean(PV[which(cycle > 800)]$R20m)
sd(PV[which(cycle > 800)]$R30m) / mean(PV[which(cycle > 800)]$R30m)
sd(PV[which(cycle > 800)]$R50m) / mean(PV[which(cycle > 800)]$R50m)
sd(PV[which(cycle > 800)]$R100m) / mean(PV[which(cycle > 800)]$R100m)

sd(PV$R5m) / mean(PV$R5m)
sd(PV$R10m) / mean(PV$R10m)
sd(PV$R20m) / mean(PV$R20m)
sd(PV$R30m) / mean(PV$R30m)
sd(PV$R50m) / mean(PV$R50m)
sd(PV$R100m) / mean(PV$R100m)

wilcox.test(PV[which(cycle > 600 & cycle <= 800)]$R30,
            PV[which(cycle > 800 & cycle <= 1000)]$R30,
            paired = FALSE, alternative = "two.sided")

g1 <- ggplot(data = PV, aes(x=cycle, y=R5m)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = R5m + R5s, ymin = R5m - R5s),
              alpha = 0.2, colour = 0, linetype = 1) + ylim(0, 1) +
  ylab('Mean Violation') + xlab('Cycle') + labs(title = "5 Replications") +
  theme(plot.title = element_text(lineheight=.8, face="bold"),
        axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

g2 <- ggplot(data = PV, aes(x=cycle, y=R10m)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = R10m + R10s, ymin = R10m - R10s),
              alpha = 0.2, colour = 0, linetype = 1) + ylim(0, 1) +
  ylab('Mean Violation') + xlab('Cycle') + labs(title = "10 Replications") +
  theme(plot.title = element_text(lineheight=.8, face="bold"),
        axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

g3 <- ggplot(data = PV, aes(x=cycle, y=R20m)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = R20m + R20s, ymin = R20m - R20s),
              alpha = 0.2, colour = 0, linetype = 1) + ylim(0, 1) +
  ylab('Mean Violation') + xlab('Cycle') + labs(title = "20 Replications") +
  theme(plot.title = element_text(lineheight=.8, face="bold"),
        axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

g4 <- ggplot(data = PV, aes(x=cycle, y=R30m)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = R30m + R30s, ymin = R30m - R30s),
              alpha = 0.2, colour = 0, linetype = 1) + ylim(0, 1) +
  ylab('Mean Violation') + xlab('Cycle') + labs(title = "30 Replications") +
  theme(plot.title = element_text(lineheight=.8, face="bold"),
        axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

g5 <- ggplot(data = PV, aes(x=cycle, y=R50m)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = R50m + R50s, ymin = R50m - R50s),
              alpha = 0.2, colour = 0, linetype = 1) + ylim(0, 1) +
  ylab('Mean Violation') + xlab('Cycle') + labs(title = "50 Replications") +
  theme(plot.title = element_text(lineheight=.8, face="bold"),
        axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

g6 <- ggplot(data = PV, aes(x=cycle, y=R100m)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = R100m + R100s, ymin = R100m - R100s),
              alpha = 0.2, colour = 0, linetype = 1) + ylim(0, 1) +
  ylab('Mean Violation') + xlab('Cycle') + labs(title = "100 Replications") +
  theme(plot.title = element_text(lineheight=.8, face="bold"),
        axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

png("exp1.png", width=1024, height=768)
multiplot(g1,g3,g5,g2,g4,g6, cols=2)
dev.off()

mean(PV[which(cycle > 800 & cycle <= 1000)]$R5)
mean(PV[which(cycle > 800 & cycle <= 1000)]$R10)
mean(PV[which(cycle > 800 & cycle <= 1000)]$R20)
mean(PV[which(cycle > 800 & cycle <= 1000)]$R30)
mean(PV[which(cycle > 800 & cycle <= 1000)]$R50)
mean(PV[which(cycle > 800 & cycle <= 1000)]$R100)


ggplot(data = mergeQS, aes(x=cycle, y=mean)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = mean + sd, ymin = mean - sd),
              alpha = 0.2, colour = 0, linetype = 1)


ggplot(data = mergeQB, aes(x=cycle, y=mean)) +
  geom_line(size = 0.5) +
  geom_ribbon(aes(ymax = mean + sd, ymin = mean - sd),
              alpha = 0.2, colour = 0, linetype = 2) +
  ylab('Mean kWh Bought') + xlab('Cycle') +
  theme(axis.title.x = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.title.y = element_text(colour = 'black', size = 16, face = 'bold'),
        axis.text.x = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.text.y = element_text(colour = 'black', size = 12, face = 'bold'),
        axis.line = element_line(colour = 'black', size = 0.5, linetype = 'solid'),
        panel.background = element_rect(fill = "transparent",colour = NA),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

limit <- 900
prop <- mean(mergePV[which(cycle >= limit)]$mean)
nv <- sum(mergeNV[which(cycle >= limit)]$mean)
nc <- sum(mergeNC[which(cycle >= limit)]$mean)
nc / (nc + nv)
nc
nv
sum(mergeQS[which(cycle >= limit)]$mean)
sum(mergeQB[which(cycle >= limit)]$mean)



# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}